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Abstract 

A general stability analysis is presented for the determination of the transition from incoherent to 
coherent behavior in an ensemble of globally coupled, heterogeneous, continuous-time dynamical 
systems. The formalism allows for the simultaneous presence of ensemble members exhibiting 
chaotic and periodic behavior, and, in a special case, yields the Kuramoto model for globally 
coupled periodic oscillators described by a phase. Numerical experiments using different types of 
ensembles of Lorenz equations with a distribution of parameters are presented. 
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I. INTRODUCTION 



Systems of many coupled dynamical units are of great interest in a wide variety of scientific 
fields including physics, chemistry and biology. In this paper, we are interested in the case 
of global coupling in which each element is coupled to all others. Beginning with the work 
of Kuramoto Jlj] and Winfree 0, there has been much research on synchrony in systems of 
globally coupled limit cycle oscillators (e.g., |3|, f|, || ]B|, |7], ||, |10[). Possible applications 
include groups of chirping cricket s|Tl]], flashing fireflies ||12|| , Josephson junction arrays [13 



semiconductor laser arrays M, and cardiac pacemaker cells []14[ . Recently, Pikovsky, et 
al. |TJJ and Sakaguchi [|nj have studied the onset of synchronization in systems of globally 
coupled chaotic systems. 

In this paper we present and apply a formal analysis of the stability of the unsynchro- 
nized state (or "incoherent state") of a general system of globally coupled heterogeneous, 
continuous-time dynamical systems. In our treatment, no a priori assumption about the 
dynamics of the individual coupled elements is made. Thus the systems can consist of el- 
ements whose natural uncoupled dynamics is chaotic or periodic, including the case where 
both types of elements are present. Our treatment is related to the marginal stability inves- 



tigation of Ref. [pu|| ; see also ||16|| . The main difference between our work and these previous 
works is that we treat an ensemble of nonidentical systems, considering both chaotic and 
limit cycle dynamics, and that our work yields growth rates as well as instability conditions. 
In addition, our treatment addresses some basic issues of the linear theory (e.g., analytic 
continuation of the dispersion function). 

The organization of the rest of this paper is as follows. The problem is formulated in 
Sec. II, and a formal solution for the dispersion relation D(s) = is given in Sec. fH\ . 
Here the quantity s governs the stability of the system (Re(s) > implies instability). 
The interpretation, analytic properties, and numerical calculation of the dispersion relation 
are discussed in Sec. |TV] along with other issues related to the treatment given in Sec. [HI. 
In Sec. [V|, we obtain D(s) for the Kuramoto model of coupled limit cycle oscillators as an 
example. Section [VT] presents illustrative numerical examples using three different ensembles 
of globally coupled Lorenz equations. In particular, these ensembles are formed of systems 
with a parameter r that is uniformly distributed in an interval [r_,r + ] with three different 
choices of r±. In the first example (Sec. |VI A| ) every element of the uncoupled ensemble 
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is periodic, but the interval [r_,r + ] includes a pitchfork bifurcation. The second example 
(Sec. |V1B| ) is for an apparently chaotic ensemble, while the third example (Sec. |V1 C| ) involves 
an ensemble that includes both chaotic and periodic elements. Finally, Sec. |V11| concludes 
the paper with further discussion and a summary of the results. 



II. FORMULATION 

We first treat the simplest case, giving generalizations later in the paper (Sec. |1V F| ). We 
consider dynamical systems of the form 

dx i (t)/dt=G(x i (t),n i ) + K(((x)). - «x(t))», (1) 

where Xj = {x\ , x\ , . . . , x\ ) T is a g-dimensional vector; G is a g-dimensional vector func- 
tion; K is a constant q x q coupling matrix; i — 1, 2, - • • , N is an index labeling components 
in the ensemble of coupled systems (in our analytical work we take the limit N — > oo, while 
in our numerical work N » 1 is finite); ((x(t))) is the instantaneous average component 
state (also referred to as the order parameter), 

«x(t)»= lim N- lJ £(Mt)), (2) 

i 

and, for each i, (xj) is the average of x« over an infinite number of initial conditions Xj(0), 
distributed according to some chosen initial distribution on the attractor of the ith uncoupled 
system 

dx i /tZt = G(x i ,n i ). (3) 

fli is a parameter vector specifying the uncoupled (K = 0) dynamics, and ((x))* is the natural 
measure [0 and i average of the state of the uncoupled system. That is, to compute ((x))*, 
we set K = 0, compute the solutions to Eq. (^), and obtain ((x))* from 

/"TO 

«x», = lim N- 1 J2l I™ r, 1 / ^(t)dt\. (4a) 

In what follows we assume that the Qi are randomly chosen from a smooth probability 
density function p(fl). Thus, an alternate means of expressing (f|a) is 

«x))„ = J xp(fi)rf/i n rffi, (4b) 
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where /in is the natural invariant measure for the system dx/dt = G(x, Q). By construction, 
((x)) = ((x))* is a solution of the globally coupled system ([[]). We call this solution the 
"incoherent state" because the coupling term cancels and the individual oscillators do not 
affect each other. The question we address is whether the incoherent state is stable. In 
particular, as a system parameter such as the coupling strength varies, the onset of instability 
of the incoherent state signals the start of coherent, synchronous behavior of the ensemble. 

III. STABILITY ANALYSIS 

To perform the stability analysis, we assume that the system is in the incoherent state, 
so that at any fixed time t, and for each i, Xj(t) is distributed according to the natural 
measure. We then perturb the orbits Xj(t) — > Xj(t) + 5xj(t), where 5xj(t) is an infinitesimal 
perturbation: 

dSxi/dt = DG(xi(t), QJSxi - K((5xi)) (5) 

where 

d 

DG(x 4 (t),a)^ = 5xi ■ — G(x i (t),fi i ). 

CXj 

Introducing the fundamental matrix Mj(t) for system (|5]), 

dMi/dt = DG • M i; (6) 
where Mj(0) = 1, we can write the solution of Eq. @ as 

5xi(t) = - f M l (t)Mr 1 (r)K((<5x)) rC /r, (7) 

J — oo 

where we use the notation ((5x)) T to signify that ((<5x)) is evaluated at time r. Note that, 
through Eq. (P), Mj depends on the unperturbed orbits Xj(t) of the uncoupled nonlinear 
system (0), which are determined by their initial conditions Xj(0) (distributed according to 
the natural measure). 

Assuming that the perturbed order parameter evolves exponentially in time (i.e., ((5x)) = 
Ae st ), Eq. © yields 

{1 + M(s)K}A = 0, (8) 
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where s is complex, and 

M(s) = ^y e-^Mi^MT 1 ^)^ . (9) 

Thus the dispersion function determining s is 

D(s) = det{l + M(s)K} = 0. (10) 

In order for Eqs. @ and (pTOj) to make sense, the right side of Eq. (|9]) must be independent 
of time. As written, it may not be clear that this is so. We now demonstrate this, and express 
M(s) in a more convenient form. To do this, we make the dependence of Mj in Eq. (||) on the 
initial condition explicit: M i (t)Mr 1 (r) = M i (t,x i (0))Mr 1 (r,x i (0)). From the definition of 
Mj, we have 

Mi(t, Xi(p))Mf\T, x,,(0)) = Mi(t - r, x,(r)) = M,(T, x,(t - T)), (11) 

where we have introduced T = t — t. Using Eq. (|TT|) in Eq. (H) we have 

M(s) = ^jf e-^MiiT^it-T^T^J . 

Note that our solution requires that the integral in the above converge. Since the growth of 
Mj with increasing T is dominated by hi, the largest Lyapunov exponent for the orbit Xj, 
we require 

Re(s) > T , T = max/ij. 

In contrast with the chaotic case where T > 0, an ensemble of periodic attractors has T = 
(for an attracting periodic orbit hi = corresponds to orbit perturbations along the flow). 
With the condition Re(s) > T, the integral converges exponentially and uniformly in the 
quantities over which we average. Thus we can interchange the integration and the average, 

POO 

M(s)= / e- sT ((M l (T,x,(t-T))))^T. (12) 
Jo 

In Eq. Ql2]) the only dependence on t is through the initial condition Xj(£— T). However, since 
the quantity within angle brackets includes not only an average over i, but also an average 
over initial conditions with respect to the natural measure of each uncoupled attractor 
i, the time invariance of the natural measure ensures that Eq. (|12|) is independent of t. 
In particular, invariance of a measure means that if an infinite cloud of initial conditions 



Xj(0) is distributed on uncoupled attractor i at t = according to its natural invariant 
measure, then the distribution of the orbits, as they evolve to any time t via the uncoupled 
dynamics (Eq. (|3])), continues to give the same distribution as at time t = 0. Hence, 
although Mj(T, Xj(t — T)) depends on t, when we average over initial conditions, the result 
(Mj(T, Xj(t — T)))* is independent of t for each i. Thus we drop the dependence of ((Mi))* 
on the initial values of the Xj and write 

POO 

M(s) = / e- sT ((M(T)»^T, (13) 
Jo 

where, for convenience we have also dropped the subscript %. Thus M is the Laplace trans- 
form of ((M))*. This result for M(s) can be analytically continued into Re(s) < T, as 
explained in Sec. [IV A| . 

Note that M(s) depends only on the solution of the linearized uncoupled system (Eq. (|6])). 
Hence the utility of the dispersion function D(s) given by Eq. ([TOD is that it determines 
the linearized dynamics of the globally coupled system in terms of those of the individual 
uncoupled systems. 

IV. DISCUSSION 

A. Analytic Continuation of M(s) 

Consider the kth column of ((M(i)))*, which we denote [((M(i)))*]fc. According to our 
definition of Mj given by Eq. (|]), we can interpret [((M(£)))*]fc as follows. Assume that for 
each of the uncoupled systems % in Eq. (0), we have a cloud of an infinite number of initial 
conditions sprinkled randomly according to the natural measure on the uncoupled attractor. 
Then, at t — 0, we apply an equal infinitesimal displacement 5k in the direction k to each 
orbit in the cloud. That is, we replace Xj(0) by Xj(0) + 5k&k, where is a unit vector in 
x-space in the direction k. Since the particle cloud is displaced from the attractor, it relaxes 
back to the attractor as time evolves. The quantity [((M))*]^ gives the time evolution of 
the i-averaged perturbation of the centroid of the cloud as it evolves back to the attractor 
and redistributes itself on the attractor. 

We now argue that ((M))* decays to zero exponentially with increasing time. We consider 
the general case where the support of the smooth density p(O) contains open regions of Q for 
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which the dynamical system @ has attracting periodic orbits as well as a positive measure 
of fl on which Eq. (|3|) has chaotic orbits. Numerical experiments on chaotic attractors 
(including structurally unstable attractors) generally show that they are strongly mixing; 
i.e., a cloud of many particles rapidly arranges itself on the attractor according to the natural 
measure. Thus, for each f2j giving a chaotic attractor, it is reasonable to assume that the 
average of Mj over initial conditions Xj(0), denoted (Mi)*, decays exponentially. For a 
periodic attractor, however, (Mi)* does not decay: a distribution of orbits along a limit 
cycle comes to the same distribution after one period, and this repeats forever. Thus, if 
the distribution on the limit cycle was noninvariant, it remains noninvariant and oscillates 
forever at the period of the periodic orbit. On the other hand, periodic orbits exist in open 
regions of fl, and, when we average over Q, there is the possibility that with increasing time 
cancellation causing decay occurs via the process of "phase mixing" [IS]. For this case we 
appeal to an example. In particular, the explicit computation of (Mi)* for a simple model 
limit cycle ensemble is given in Sec. 0. The result is 



(Mi), 



cos flit — sin flit 
sin Qit cos flit 



and indeed this oscillates and does not decay to zero. However, if we average over the 
oscillator distribution p(fl) we obtain 



«M», 



c(t) -s(t) 
s(t) c(t) 



where c(t) = J p(Q)cosQtdQ and s(t) = J p(fl) sin QtdQ. For any analytic p(fl) these 
integrals decay exponentially with time. Thus, based on these considerations of chaotic 
and periodic attractors, we see that for sufficiently smooth p(fl), there is reason to believe 
that ((M))*, the average of Mi over Xj(0) and over decays exponentially to zero with 
increasing time. Conjecturing this decay to be exponential, | ((M(i)))*| < kc~& for positive 
constants k and £, we see that the integral in Eq. (^) converges for Re(s) > — £. This 
conjecture is supported by our numerical results in Sec. [VT|. Thus, while Eq. (|T3]) was 
derived under the assumption Re(s) > T > 0, using analytic continuation, we can regard 
Eq. (|T3|) as valid for Re(s) > — £. Note that, for our purposes, it suffices to require only that 
||((M(i)))*|| be bounded, rather than that it decay exponentially. Boundedness corresponds 
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to £ = 0, which is enough for us, since, as soon as instability occurs, the relevant root of 
D(s) has Re(s) > 0. 



B. Numerically Approximating M(s) by use of Eq. (^) 

We can envision the following numerical method for finding ((M(£)))*. First approximate 
the natural measure of each attractor i by a large finite number of orbits initially distributed 
according to the natural measure. For each initial condition, obtain Xj(t) from Eq. (|3]). Then 
use these solutions in DG and solve Eq. @. Finally, average the resulting matrix solutions 
Mj over the orbits. While this may look attractive, it can present difficulties in the case 
of chaotic orbits. In particular, chaos implies that the individual Mj diverge exponentially, 
while the average over the natural invariant measure (Mj)* decays. That is, when we average 
over the natural invariant measure, the exponential divergence of the individual Mj(i) cancel 
to yield decay. Numerically, however, we average over a large but finite number of orbits. 
For early time, one can expect that this will give a good approximation to (Mj)*. However, 
as time goes on, the decay of (Mj)* implies that the cancellation must become more and 
more accurate because the individual Mj become larger and larger. Thus, eventually, any 
numerical approximation using a finite number of orbits must diverge. The question is, 
can one obtain results by this method that are accurate for long enough time to provide a 
useful basis for approximating M(s). We expand and illustrate this issue in greater detail 



in Sec. VI 



A variant of the above numerical technique is to obtain ((M))* by working directly with 
the uncoupled nonlinear equations (0). We use a large finite number of initial conditions 
chosen randomly with respect to p(fl) and the natural invariant measure. These initial 
conditions are then all displaced, Xj(0) — > Xj(0) + where A& = A^a^. and A& is small. 
Denoting the solutions of Eq. (|]) with these displaced initial conditions x^(t, A&), we then 
approximate the quantity (((x))* — ((x'(i, A&)))), which represents the relaxation of the 
measure's centroid after a displacement A&. In the limit A& — > 0, we have that the kth 
column of ((M(t)))* is 

[((M(t)))*] fe = A fc 1 (((x))*-((x'))), (14) 
and Eq. ( |14|) is used, with small A&, as an approximation. Again, practical numerical issues 
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exist for this technique. In particular, must be small for linearity to be approximated, 
but this makes the cancellation of ((x))* — ((x')) stronger, which, in turn, necessitates using 
many initial conditions. Also, as time increases, ((M(£)))* decreases, and fluctuations of 
((x))* — ((x')) due to the finite number of initial conditions can overwhelm the computation 
of the coherent relaxation to the attractor (see Sec. |VJ). 

C. Numerical approximation of M(s) as the Response to exp(si)l 

Since ((M))* is the response to an impulse drive, its Laplace transform multiplied by e st , 
M(s)e si , is the response to a drive e st t. We now show this more formally. First we note that 
Eq. (|6|) with the initial condition M, = 1 at t — can be written in the impulse response 
form, 

dMi/dt = BG ■ Mi + 5(t)l, 

where 5(t) is a delta function, and Mj satisfies the initial condition Mj = at t = — oo. 
Shifting time by £q, we have 

-Mi(t - t , Xi(t )) = DG • Mi(t - t , Xi(* )) + 6(t -t )l, (15) 

where we have explicitly included the dependence of Mj on time and on the initial state 
Xj(to) of the unperturbed orbit Xj(t). Multiplying by e st °dt Q and integrating over all to we 
obtain 

dM.i/dt = DG • Mj + e st l, (16) 

where 

Mi= f e s ' Mj(t-t ,x,(io))^o, (17) 

J — oo 

which converges at the lower limit provided Re(s) is sufficiently large. Introducing T = t— 1 , 
averaging, and, as before, interchanging the order of the average and the integral, we have 
that the response to e st t is, as claimed, ((M))* = e st M(s). This suggests the following 
numerical technique for finding M(s). Solve 

—tt- = G(xf' s) ,fij) + aJ (18) 
ai L e sin uit, 
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where s = o — iu, A^ = Ak&k, and A^ is real. For large t, but A^e "* still small throughout 
the interval (0,t), we can regard the average response as approximately linear. Thus the 
kth column of M(s) is 

[M( S )] fc ^A fc -V st [((x)),-((x))], (19) 

where Xj = x- — iic^ . Numerically ((x)) can be approximated using a large finite number 
of orbits. In Ref. |16|| , a technique equivalent to this with s taken to be imaginary (s = 
— iu) is used to obtain the marginal stability condition (see also Ref. |17|]). In Sec. [Vj we 



compare the numerical efficacy of this technique and of the techniques discussed in Sec. IV B 



20fl . The reasoning in Ref. |L6| is heuristic, and, adapted to our setting | 2lR , it goes as 



follows. Numerically, it is observed that as the coupling is varied, a Hopf bifurcation occurs. 
Thus, for conditions just past the bifurcation, the order parameter variation is sinusoidal, 
((x))* — ((x)) ~ e~ luJt . Using this as the drive in the nonlinear equation, as in Eq. ([18]), and 
computing M(— iu) as above, self-consistency then yields {1 + M(— iu) ■ K}A = 0. For the 
case of a coupling matrix with one nonzero element located on the diagonal [i.e., K\\ = k 
and Kij = if ^ (1,1)], the consistency condition then gives 1 + M\i{— iui)k = 0. 
Setting the real and imaginary parts of this equation equal to zero determines the value of 



the frequency at instability onset and the critical value of the coupling constant k p2 



D. The Distribution Function Approach 

Much previous work has treated the Kuramoto problem and its various generalizations 
using a kinetic equation approach^, f|, [|, H, H, H HH, [TTj, We have also obtained 



our main result, Eq. (pT|) for D(s), by this more traditional method. We briefly outline the 
procedure below. 

Let -F(x, Q, t) be the distribution function (actually a generalized function) such that 
F(x, Q, t)dxdQ is the fraction of oscillators at time t whose state and parameter vectors lie 
in the infinitesimal volume d'x.dVL centered at (x, Q). Note that J Frfx is time independent, 
since it is equal to the distribution function p(fl) of the oscillator parameter vector. The 
time evolution of F is simply obtained from the conservation of probability following the 
system evolution, 

BF B 

— + — • [(G(x, O) + K ■ («x», - «x»))F] = 0, (20) 
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where 



(21) 




((x))* = / / dxrffixFo, (22) 



and Fq = -Fo(x, Q) = /(x, Q)p(Q), in which /(x, Q) is the density corresponding to the 
natural invariant measure of the uncoupled attractor whose parameter vector is Q. Thus 
/(x, Q), which is a generalized function, formally satisfies 

|--[G(x,Q)/(x,O)] = 0. (23) 



Hence, F = F is a time-independent solution of Eq. (|20|) (the "incoherent solution"). We 
examine the stability of the incoherent solution by linearly perturbing F, F = Fq + 5F, to 
obtain 

^_ + _ . [ G (x, n)5F - K((5x))F ] = (24) 



«5x)) = / / dxdnxSF. (25) 




We can then introduce the Laplace transform, solve the transformed version of Eq. (p4|). 
and substitute into Eq. ( p5|) to obtain the same dispersion function D(s) as in Sec. |HI[ The 
calculation is somewhat lengthy, involving the formal solution of Eq. (^4|) by integration 
along the orbits of the uncoupled system. We will not present the detailed steps here, since 
the result is the same as that derived in Sec. JTJ, where it is obtained in what we believe is 
a more direct manner. 

The distribution function approach outlined above is similar to the marginal stability 
treatment of Ref. |L7| for identical globally chaotic maps. In that case s — > —iu, the 



Frobenius- Perron equation plays the role of Eq. (^0|), and the average over parameters is not 
present. 

We note that the computation outlined above is formal in that we treat the distribution 
functions as if they were ordinary, as opposed to generalized, functions. In this regard, 
we note that /(x, Q) is often extremely singular both in its dependence on x (because 
the measure on a chaotic attractor is typically a multifractal) and on Q (because chaotic 
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attractors are often structurally unstable). We believe that both these sources of singularity 
are sufficiently mitigated by the regularizing effect of the averaging process over (x, fl), and 
that our stability results of Sec. [TTT] are still valid. This remains a problem for future study. 
We note, however, that for structurally unstable attractors, a smooth distribution of system 
parameters p(fi) is likely to be much less problematic than the case of identical ensemble 
components^, |T7j ], p(&) = S(Q — O). In the case of identical structurally unstable chaotic 
components, an arbitrarily small change of Cl can change the character of the base state 
whose stability is being examined. In contrast, a small change of a smooth distribution p(Q) 
results in a small change in the weighting of the ensemble members, but would seem not to 



cause any qualitative change. We remark that the numerical test cases we study in Sec. [VI 
are all structurally unstable. Nevertheless, they all agree well with the theory. 

E. Bifurcations 

It is natural to ask what happens as a parameter of the system passes from values cor- 
responding to stability to values corresponding to instability. Noting that the incoherent 
state represents a time independent solution of Eq. flU), we can seek intuition from stan- 
dard results on the generic bifurcations of a fixed point of a system of ordinary differential 



equations (||23||; see also |17|). There are two linear means by which such a fixed point can 
become unstable: (i) a real solution of D(s) = can pass from negative to positive s values, 
and (ii) two complex conjugate solutions, s and s*, can cross the imaginary s-axis, moving 
from Re(s) < to Re(s) > 0. 

In reference to case (i), we note that the incoherent steady state always exists for our 
formulation in Sec. II. In this situation, in the absence of a system symmetry, the generic 
bifurcation of the system is a transcritical bifurcation (Fig. 0(a)). In the presence of 
symmetry, the existence of a fixed point solution with ((x))* — ((x)) nonzero may imply the 
simultaneous existence of a second fixed point solution with ((x))* — ((x)) nonzero, where 
these solutions map to each other under the symmetry transformation of the system. In this 
case the transcritical bifurcation is ruled out, and the generic bifurcation is the pitchfork 
bifurcation, which can be either subcritical (Fig. 0(b)) or supercritical (Fig. |l](c)) . 

In case (ii), where two complex conjugate solutions cross the Im(s) axis, the generic 
bifurcations are the subcritical and supercritical Hopf bifurcations. (In this case we note 
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that although the individual oscillators may be behaving chaotically, their average coherent 
behavior is periodic.) 

In our numerical experiments in Sec. |VT| we find cases of apparent subcritical and su- 
percritical Hopf bifurcations, as well of what we believe is a subcritical pitchfork 
bifurcation. The reason we believe it is a pitchfork rather than a transcritical bifurcation 
is that our globally coupled system is a collection of coupled Lorenz equations. Since the 
Lorenz equations 

dx^/dt = a(x^-x^) 

dxW/dt = rx« - x& - X W X (3) (26) 
dx^/dt = -6x (3) + x (1) x (2) 

have the symmetry (x^\ x^^x^) — ► (— x^ l \ —x^ 2 \ x^), and since the form of the coupling 
used in Sec. [VT| respects this symmetry, the transcritical bifurcation is ruled out. 

F. Generalizations 

One generalization is to consider a general nonlinear form of the coupling such that we 
replace system flip by 

dxi/dt = G(xi,J2i,y) 
y = «x»„-«x» 

and the role of the uncoupled system (analogous to Eq. (^)) is played by the equation 

dXj/dt = G(x i ,n i ,0). (28) 
In this more general setting, following the steps of Sec. |IIJ yields 

D(s) = det{l + Q(s)}, (29) 

where 

poo 

Q0)= / dTe- st ((M(T)D y G(x,ft,0)>>*. 
Jo 

A still more general form of the coupling is 

dxi/tft=G(x i ,fi i ,((x))). (30) 
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For Eqs. (|27D and (|l|), a unique incoherent solution ((x))* always exists and follows from 
Eq. (^) by solving the nonlinear equations for each x,(0) with y = (((x))* — ((x))) set 
equal to zero. In the case of Eq. (|30]), the existence of a unique incoherent state is not 
assured. By definition, ((x)) is time independent in an incoherent state. Thus replacing 
((x)) in Eq. ( |50| ) by a constant vector u, imagine that we solve Eq. fl5UD for an infinite 
number of initial conditions distributed for each i on the natural invariant measure of the 
system, chti/dt = G(xj,f2j,u), and then use Eq. (§) to obtain the average ((x)) u . This 
average depends on u, so that ((x)) u = F(u). We then define an incoherent solution ((x))* 
for Eq. (|30|) by setting ((x)) u = u = ((x))*, so that ((x))* is the solution of the nonlinear 
equation 

«x». = F(«x».). 

Generically, such a nonlinear equation may have multiple solutions or no solution. In this 
setting, if a stable solution of this equation exists for some paramter k < k c , then the 
solution of the nonlinear system (|30| ) (with appropriate initial conditions) will approach it 
for large t. If now, as k approaches k c from below, a real eigenvalue approaches zero, then 
k = k c generically corresponds to a saddle-node bifurcation. That is, an unstable incoherent 
solution merges with the stable incoherent solution, and, for k > k c , neither exist. In this 
case, loss of stability by the Hopf bifurcation is, of course, still generic, and the incoherent 
solution continues to exist before and after the Hopf bifurcation. D(s) for Eq. (P0|) is given by 
Eq. (|29|) with D y G replaced by — D(( X ))G evaluated at the incoherent state (((x)) = ((x))*) 
whose stability is being investigated. 

Another interesting case is when the coupling is delayed by some linear deterministic 
process. That is, the ith oscillator does not sense ((x)) immediately, but rather responds 
to the time history of ((x)). Thus, using Eq. fl27[) as an example, the coupling term y is 
replaced by a convolution, 

y(f)= f rft'A(t-t')-(«x»*-«x» t ,). 

J — oo 

In this case a simple analysis shows that Eq. (|29|) is replaced by 

D(s) = det{l + Q(s)-A(s)} 

where 

POO 

A(s) = / dte- st A(t'). 
Jo 
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The simplest form of this would be a discrete delay 

in which case A(s) = te~~ vs . (The case of time delayed interaction has been studied for 
coupled limit cycle oscillators in Refs. || [7], |J.) 

In addition to these generalizations, others are also of interest. For example, the inclusion 
of noise and coupling "inertia" is studied in the limit cycle case in Ref . || . 



V. THE KURAMOTO PROBLEM 



As an example, we now consider a case that reduces to the well-studied Kuramoto prob- 
lem. We consider the ensemble members to be two dimensional, Xj = (xi(t),yi(t)) T , and 
characterized by a scalar parameter fli. For the coupling matrix K we choose fcl. Thus 
Eq. (|1|) becomes dxi/dt = G {x) {xi, y h fli) + &((<»)* - ((a:))), dyt/dt = G^^x^y^Qi) + 
k({{y))* ~ ((y)))- We assume that in polar coordinates (x = rcos9,y = rsin^), the uncou- 
pled (k — 0) dynamical system is given by 

dOi/dt = n i} (31) 



dri/dt = (r - r^/r, 



(32) 



where f^r <C 1. That is, the attractor is the circle = r , and it attracts orbits on a 
time scale r that is very short compared to the limit cycle period. For flir <d 1 it will 
suffice to calculate Mj(t) for t ^> r. To do this, as shown in Fig. |2|, we consider an initial 
infinitesimal orbit displacement A oi = a. x dx i + dL y dy oi where a. XtV are unit vectors. In a 
short time this displacement relaxes back to the circle, so that for (2-jr/fl) > t > r we have 
r = r , 9 = 9 oi , Aj-(t) ~ A+a e , where oi is the initial value 0j(O), a g is evaluated at 6^(0), 
and A+ = — sin 9 oi dx oi + cos 9 oi dy oi . For later time t ^> r, we have r = r , ^(t) = 9 oi + flit 
and Aj(t) = A^a#, with a e evaluated at 6i(t). In rectangular coordinates this is 

" sm(9 oi + flit) sin 9 io - sm.(9 oi + flit) cos io 1 |~ dx oi 
. — cos(9 oi + fltf) sin O j cos(0 O j + flit) cos O j J L dyoi 



dxi(t) 




. dyi(t) _ 
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By definition, the above matrix is Mj appearing in Sec. pTIl Averaging Eq. (0) over the 
invariant measure on the attractor of Eqs. fl3~T|) and (p2|) implies averaging over O i. This 
yields 



(Mi) f 



cos flit — sin flit 
sin cos flit 



(33) 



Averaging the rotation frequencies fli over the distribution function p(fl) and taking the 
Laplace transform gives M(s), 



M(s) 



(q + + q_) i(q + - q_ 
-i(q + - qJ) (q + + q„) 



(34) 



where 



p(fl)dfl 



(35) 



4\sTifl/ Q 4: J -co s =F ^ 
and, in doing the Laplace transform, we have neglected the contribution to the Laplace 
integral from the short time interval < t < 0(r) (this contribution approaches zero as 
fir -> 0). Using Eqs. and (|35|) in Eq. (Jl0|) then gives D(s) = D + (s)D_{s), where D±(s) 
is the well-known result for the Kuramoto model (e.g., [|HJ), 



£>±(s) = l + 



A- 



s ± ifl 



, ite(s) > 0, 



(36) 



and D±(s) for Re(s) < is obtained by analytic continuation P, |T^]. Note that the property 
D±(s) = D T (s^), where f denotes complex conjugation, insures that complex roots of D(s) = 



Djs)D(s) 



come in conjugate pairs. 



VI. NUMERICAL EXPERIMENTS 



In this section, we illustrate and test our theoretical results using three different ensembles 
of globally coupled heterogeneous Lorenz oscillators. The Lorenz equations are given in 
Eq. ( p6|) . For our numerical experiments, we set a = 10 and b = 8/3, and the ensembles 
are formed of systems with the parameter r uniformly distributed in an interval [r_ , r + ] . 
We consider three different cases: an ensemble of periodic oscillators containing a pitchfork 
bifurcation (r_ = 150 and r + = 165), an apparently chaotic ensemble (r_ = 28 and r + = 52), 
and an ensemble with mixed chaotic and periodic oscillators (r_ = 167 and r + = 202). 

16 



As previously stated, the dispersion function D(s) given by Eq. ( [L0|) depends only on the 
solution of the linearized uncoupled system, and D(s) in turn determines the linear stability 
of the incoherent state of the globally coupled system. To demonstrate this numerically, we 
consider the simple case in which the coupling matrix K has only one nonzero diagonal 
element, i.e., K\\ = k and Kij = for all 7^ (1,1)- For each type of ensemble, 

a large number (N > 10 4 ) of Lorenz equations (Eq. (PSD) were initialized with random 
initial conditions chosen within their respective basins of attraction. The Lorenz equations 
were integrated using the standard 4th-order fixed-time-step Runge-Kutta method. Each 
element in the ensemble was independently integrated for a sufficiently long but random 
time to ensure that the ensemble was essentially on the attractor and was not initiated 
in a coherent state. Since the number of oscillators iV is large, we choose a simpler form, 
((x^)} = N~ x xf \ for the order parameter defined in Eq. (fj). With N sufficiently large, 
the average (a^ ) over the natural measure for a given system i can be absorbed into the 
summation over i. 

In the numerical experiments below, we will consider the following time averaged quantity 

as a measure of the coherence of the order parameter: 

r r t+T -1 1/2 

x T = T" 1 J ((x {1 \t'))) 2 dt' , (37) 

where t is sufficiently long so that the ensemble has achieved its time asymptotic behavior, 
and T is sufficiently long so that xt is essentially independent of T. Note that the symmetry 
of the Lorenz equations under (x^\ x^ 2 \ x^) — > (-x' 1 ',-^ 2 ',^) implies that ((x^))* = 
for N = 00 when the initial conditions are distributed symmetrically in space. Near 
k — 0, the ensemble is weakly coupled and the incoherent solution is expected, i.e., xt ~ 0. 
(Although ((x^)) is time independent in the infinite iV limit for the incoherent state, ((x^)) 
will exhibit small fluctuations in time for finite N.) As the magnitude of k increases, 
transitions to different coherent states where xt > were observed (see Figs. p],|TT], and 
|18| below). With the three ensembles we have chosen, we have observed super-critical Hopf 
(the right transition in Fig. |5| and the left transition in Fig. [18]), sub-critical Hopf (the left 
transition in Fig. |5| and the right transition in Fig. 18), and a subcritical pitchfork (Fig. [TTJ) 
bifurcation from the incoherent state. 

Theoretical predictions for the critical coupling strength k* for each of these bifurca- 
tions can be obtained by estimating M n (s) from the corresponding uncoupled ensemble. In 
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particular, we consider the marginal stability condition described in Sec. [IV C. First, we 



numerically integrate Eq. ( |T8| ) with a = (s = —iu) and A = (A, 0,0). With A chosen to 
be sufficiently small, the average response is linear and the (1, 1) element of M(— iuj) is 

Mu(-iw) « A- 1 e-- t [((x«)), - (38) 

where ((x^)) is defined in Eq. ([19]) and can be obtained numerically from an ensemble of 
orbits solved from Eq. (|1|). Since we are using a finite number of elements in the ensemble, 
there will be statistical noise in the ensemble averages calculated on the right hand side of 
Eq. (|38D; this can be minimized by averaging over time, i.e., 

M„(-^) « (TA)- 1 f ' e-*"[((xW)), - ((x W )))dt. (39) 

Jo 

With the coupling matrix K being nonzero only in its (1,1) position, Eq. (|10D yields 



l + M n (-iu)k = 0. (40) 

The real and imaginary parts of Eq. ( |4"0"{ ) provide two equations that can be used to determine 
both the critical value of the coupling constant k* and the frequency uj* at the onset of in- 
stability. In particular, the imaginary part of the equation, Re[Mn(— iuj)] = 0, can be solved 
for uj* (note that there may be multiple roots). The real part then yields the corresponding 
critical coupling k* = — [Mu(- za;*)] -1 . To determine which of the possibly multiple roots 
for uj* are relevant, we note that as k increases or decreases from zero, a critical value k* is 
encountered at which the incoherent state first becomes unstable. Hence we are interested 
in obtaining the smallest values of \k*\ for both positive and negative k*. (For clarity, we 
will denote the negative critical value by — \k*\.) Accordingly, the relevant uj roots are those 
yielding the largest value of |Mn(— for both Mu(-iuj) > and Mu(—iu) < 0. Denot- 
ing the corresponding uj roots by uj* and uj£, respectively, it is expected that the incoherent 
state is stable in the range —\k*\ < k < k%, where 

k* afi = -l/M n H^), -|*:| < < k* b , (41) 

and that, as k increases through kl (or decreases through — |/c*|), instability ensues. 

Growth rates and frequency shifts from uj* can also be simply obtained for k near k*. 
Setting k = k* + 5k and s = —i{uj* + 5uj) + 7 in the dispersion relation 1 + kM\\{s) = 0, 
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and expanding for small 5k, 5u and 7, we obtain 

_ 5k dIm[M u (-iu)]/du 
7 ~~(fc*) 2 {dMni-t^/dul 2 ( } 

(and a similar equation for 5u), where the expression on the right side of Eq. (|42| ) is to be 
evaluated at uo = to*. Thus, instability implies that <9ReMn(— iuo*)/duo < if k*, 5k > and 
dReM n (-iu*)/duj > if k* , 5k < 0. 

These growth rates can be estimated numerically by the following procedure. For a given 
ensemble, we initiate the system in its incoherent state by setting the coupling to zero and 
integrating for a sufficiently long time. The coupling is then switched on to a value less than 
— |A; a |* or larger than kl, where the incoherent state is unstable and the order parameter 
({x^ 1 ')) begins to grow exponentially. If the transition is a pitchfork bifurcation, we expect 
((x^))(t) ~ e 7 * for ((x^)) sufficiently close to the incoherent state. The growth rate 7 
can be obtained by measuring the slope of the graph ln((x^^))(t) vs. t. If the transition 
is a Hopf bifurcation, the growth of the order parameter will be modulated by a sinusoidal 
function. In this case, the envelope of the oscillating order parameter grows exponentially 
and the growth rate can be extracted by measuring the slope of the logarithm of this envelope 
function versus time. 

A. Periodic Ensemble 

We first consider an ensemble of Lorenz oscillators with r_ = 150 and r + = 165. In 
this range of parameters, the Lorenz equations yield stable periodic orbits. As r decreases 
through a critical value r c « 154.7, the system goes through a pitchfork bifurcation in which 
a single periodic orbit symmetric under (x^\ x^ 2 \ x^) — > (— x^\ —x^ 2 \x^) bifurcates into 
a pair of asymmetric periodic orbits. The range of dynamics for the (uncoupled, k = 0) 
Lorenz equation in this parameter range is illustrated in the bifurcation diagram in Fig. |3|. 
This bifurcation diagram is constructed by plotting the maxima of the function x^ 3 '(t) in 
the Lorenz equation for t sufficiently large so that any transient is minimized. To further 
illustrate the pitchfork bifurcation, phase space plots of the limit cycle attractors at r = 165 
and r = 150 are shown in Figs. [| (a) and (b), respectively. 

Figures || (a) and (b) show plots of xt as a function of the coupling strength k. 60, 000 
oscillators were used. Data are plotted corresponding to the cases in which k decreases 
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(black squares) and k increases (grey circles). As expected, xt is practically zero (to order 
0(N~ 1 ^ 2 )) for k near 0. As k increases through kl = 0.95 ± 0.1, the incoherent state 
destabilizes and xt increases from zero. Similarly, as k decreases through —\k* a \ = —0.70 ± 
0.04, the incoherent state again destabilizes. The transition at — \k*\ is hysteretic, while the 
transition at kl is not. It is also beneficial to examine the time evolution of the instantaneous 
order parameter ((x(t))) near the onset of these transitions. This is shown in Figs. ^| (a) 
and (b), in which ({x^(t)}} versus ((x^(t))) is shown before (grey) and just after (black) 
the transitions at (a) —\k* a \ and (b) k%. These transitions are apparently subcritical and 
supercritical Hopf bifurcations, respectively. The spread in the trajectories is due to the 
finiteness of N; we find that decreasing N increases the spread. In the following, we will 
investigate the oscillation frequency and growth rate of the order parameter near these 
transition points. 

Using the frequency response method described at the beginning of this section and 
in Sec. |IVC| , we calculated Mu(—iu) as a function of uj using an ensemble of uncoupled 
elements. We plot both the real (black) and imaginary (grey) parts of Mu(-iw) in Fig. [7] 
(a). (For comparison, Fig. [7] (b) shows the results of the linear displacement method of 
Sec. [IV jj| ; see the discussion below.) For these curves, we used a forcing strength A = 0.05 
and N = 20,000 in our calculations. As one can see, Re[Mn(— iuj)\ crosses zero more 
than once, and each root corresponds to a possible solution for uj*. Note that the maxima of 
Re[Mn(— iu)] are attained very near, but not necessarily at, these uj* roots. The two critical 
values —\k'l\ and kl are predicted by Eq. fl4~i~D with u* b corresponding to the largest values 
of |Re[Mn(— ioj*)]\ for which Re[Mn(— = 0. These values are indicated by the dotted 
lines in the figure, and yield predictions of — \k* a \ = —0.72 ±0.05 and kl = 0.93 ±0.03. These 
predictions agree well with the critical transitions observed in our numerical experiments 
using the full nonlinear system (see above). In addition, the predicted frequency at the 
supercritical bifurcation at kl is col ~ 21.4. Fig. || shows the power spectrum of ((x^(t))) 
for k = 1.2, i.e., slightly greater than kl; this figure reveals a prominent peak at a frequency 
of approximately 21.4, in agreement with the predicted value of a£. 

Since the elements within this ensemble are not chaotic, the Mj(t) do not diverge in 
time, and we expect the linear displacement method described in Sec. [IV B| for estimating 
((M(t)}}# to work well. The ensemble average ((M(i)})* should decay due to "phase mixing," 
as in the Kuramoto example (see Sees. [IV A| and[V|). A graph of ((Mn(i))}* for the periodic 
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ensemble is plotted in Fig. || As one can see, ((Mn (£)))* exhibits complex oscillatory 
behavior as it decays to zero, where small fluctuations presumably due to finite iV are 
evident. To obtain Mq(— za;), we set s = —iu in the Laplace transform of ((Mil (£)))*. 
The real (black) and imaginary (grey) parts of M n (— iu;) (black) obtained by this method 
are plotted in Fig. [7] (b). These graphs generally agree with the graphs obtained using the 
frequency response method (shown in Fig. |7| (a)) except near the roots of Re[Mn(— ioJ*)] = 0, 
where the peaks were not well resolved. Attempts to improve the frequency resolution of the 
Laplace transform requires a calculation of ((Mn (£)))* for longer time. However, fluctuations 
due to the finite number of ensemble elements prevent the accurate calculation of the decay of 
((Mm (£)))* to zero for large times. Thus, N must be increased, and practical considerations 
limit the usefulness of this method (although we note that for this example, the method 
does yield good values for u* and u£). 

Similarly, an accurate measurement of the growth rate of the mean field requires very 
large ensembles and extremely long transients due to the weak phase mixing, and again 
we found this calculation to be numerically impractical. Thus, we demonstrate our growth 
rate calculations only in the computationally more feasible cases considered below, i.e. the 
chaotic and mixed ensembles. 



B. A Chaotic Ensemble 



We now consider an ensemble of Lorenz equations with r_ = 28 and r + = 52. From the 
bifurcation diagram (see Fig. |H]), the ensemble seems to consist of predominantly chaotic 
systems |[24|| . Within this range of parameters, the positive Lyapunov exponent varies 



between approximately 0.904 and 1.323. 

Once again, we examined the destabilization of the ensemble's incoherent state by plotting 



xt as a function of k. One can see in Fig. O that this chaotic ensemble has a hysteretic 



transition at —\k*\ = — 5.56±0.01. On the positive side, the incoherent state is stable up to 
the largest k value tested (k = 7). Examining the temporal dependence of the instantaneous 
order parameter ((x(t))) near the transition at — \k*\, we find that the order parameter jumps 
to one of two stable fixed points on opposite lobes of the Lorenz attractor (see Fig. |l2|) . As 
we have discussed previously (Sec. [TVE|) , this subcritical transition is expected to be a 



pitchfork bifurcation rather than a transcritical bifurcation due to the intrinsic symmetry 
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of the Lorenz equations. 

We calculated Mn(— iu) as a function of to by examining the uncoupled ensemble under 
periodic perturbation. For this case, we chose the forcing strength A to be 2 and N = 20, 000. 
(We varied the value of A by an order of magnitude from 0.5 to 5 and the result does not seem 
to change significantly; this indicates that the perturbation is sufficiently linear.) Figure 
shows the real and the imaginary parts of Mn(-iuj) versus u for this case. As compared 
with the previous example, the frequency response curve is simpler. Re[Mn(— iu)] crosses 
zero only at uj* = 0, where Re[Mn(— iw)] has a prominent peak. Using Eq. fl41[), this gives 
a critical coupling value of —\k*\ = —5.57 ± 0.15. This result agrees well with the threshold 
of instability for the incoherent state observed in the globally coupled ensemble. 

We have also attempted to calculate Mu(-iw) for this (chaotic) ensemble using the other 
two methods described in Sec. IV B| . These are: the linear method, in which the linearized 
equation [Eq. @] is solved for Mu(t) and the result is averaged, and the impulse-response 
method, in which the orbits on the attractor are displaced by a small amount in the i' 1 ' 
direction and the rate of decay back to zero is measured. The results from these methods 



are included in Fig. 13 with filled and open circles, respectively. While all methods agree 
reasonably well for uj > 2.5, the important narrow peak at uj = is missing from the results 
of both the linear and the impulse-response methods. 

This discrepancy can be understood by observing that the peak at uj = represents 
long-time dynamics. In particular, the half- width of this peak has Auj ~ 1, corresponding 
to a decay time scale of 1/ Auj m 1. In contrast, the spectrum, with this peak deleted, has 
a half-width of Auj « 8, corresponding to a much shorter time scale of approximately 0.1. 
The linear and impulse-response methods apparently resolve the short time scale well, but 
fail to adequately resolve the longer time scale. This is due to the problem that we have 
discussed in Sec. [IV B| . For the linear method, the individual Mj(£) grow exponentially in 
time, and hence the ensemble average ((M(t))} # requires a delicate canceling in order to 
remain valid for large time. Figure [14] shows a graph of ((Mn(t))}* for the linear method 
in grey. ((Mn (£)))* initially decays exponentially, but for t > 0.7, it begins to grow as the 
balanced canceling breaks down due to the fact that only a finite number of elements is used 
in the calculation. Thus, when obtaining the Laplace transform, we only integrated over the 
reliable range, i.e. < t < 0.7. This had the effect of leaving out the slower decay, which is 
vital for determining the critical coupling strength for the onset of instability in this case. 
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In contrast, when ((Mn(i)))* is measured using the impulse response method, it does not 
ultimately diverge exponentially. However, its exponential decay is masked by fluctuations 
for t > 0.7, again due to finite N; see the black curve in Fig. 

We found that the frequency response method is more reliable because the temporal 
averaging effectively reduces statistical noise. Therefore, we were able to obtain a good esti- 
mate of Mn(—iu) with only a moderate number of oscillators. The cost for these improved 
statistics is that each calculation is for only one particular value of u. This is in contrast to 
the impulse response method, in which the Laplace transform of ((Mn (£)))* gives the entire 
dependence of Mn(— iuf) on uo at once. Some of the comparative advantages and drawbacks 
among the three numerical methods in calculating Mu(-iuj) can be clearly seen in this 
example. 

The growth rate of the incoherent state, when it first becomes unstable, can be estimated 
from dIm[Mxi(—iu)]/du at u = uo* = 0. According to Eq. fl42|) , this growth rate (7) is a 
linear function of 5k for k close to k* . To verify this, we obtained growth rates for various 
values of k by first initializing the ensemble in the incoherent state, and then fitting a line 
to the graph of the ln((x^)) versus time. Since the transition is subcritical, only the initial 



growth rate is measured. Figure 15 shows the typical behavior for k slightly beyond the 



critical value. A plot of estimated growth rates versus 5k is shown in Fig. [TBI. The predicted 



slope, calculated from the frequency response method using Eq. (|TH), is —0.29 ± 0.02. This 
agrees well with the measured growth rates. 

C. A Mixed Chaotic Ensemble with Periodic Windows 

For our last example, we consider an ensemble which contains both chaotic and periodic 
oscillators (r_ = 167 and r + = 202). As one can see from the bifurcation diagram in Fig. [17], 
there is a prominent period-four window near r = 181. Thus, the chaotic attractors in this 



ensemble are expected to be structurally unstable. Figure 18 shows a plot of Xt clS 8b function 
of k. The incoherent state becomes unstable as k increases through kl = 1.75 ± 0.05 and 
as k decreases through —\k*\ = —0.68 ± 0.03. For the left (negative) transition, the order 
parameter ((x(t))) becomes periodic, and the amplitude of its oscillation gradually increases 
as k moves beyond its critical value at — \k*\. Thus, this transition is a supercritical Hopf 
bifurcation. The frequency of oscillation of the order parameter near this transition was 
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estimated from the power spectrum of ((x^(t))) to be u* ~ 12.2. The transition at k^ 
appears to be a (hysteretic) subcritical Hopf bifurcation. Phase portraits for k on either 
side of the two transitions are shown in Figs. [19] (a) and (b). 

As before, these two transitions can be predicted from Mu(-iuj) calculated from the 
uncoupled ensemble. Figure |2(] is a graph of the real and the imaginary parts of Mn(-iuj) 
obtained using the frequency response method (A = 0.7 and iV = 20,000). As in the 
periodic ensemble case, the maxima of Re[Mn(— iu)] occur very near, but not exactly at 
the uj roots of Re[Mn(— io;)] = 0. Using Eq. fl4"ip, the values of Re[Mn(— iw* b )] near 
the two biggest peaks give -\k*\ = -0.72 ± 0.09 and k% = 1.64 ± 0.08. The predicted 
transition frequency associated with the supercritical transition at —\k*\ is approximately 
12.2. These predictions agree well with the observed quantities obtained using the fully 
nonlinear, globally coupled ensemble. 

We have also compared the actual growth rate obtained from the globally coupled ensem- 
ble with its predicted value calculated from Mn(—iu) using the same procedure described 
above. Figure |2l| is a graph of 7 vs. 5k for the transition at — \k*\. The predicted slope, 



calculated using the frequency response method using Eq. (|18|), is —0.26 ± 0.05; this agrees 
well with the measured growth rates. 



VII. CONCLUSION 



We have presented a general formulation for the determination of the stability of the 
incoherent state of a globally coupled system of continuous time dynamical systems. This 
formulation gives the dispersion function in terms of a matrix M(s) which specifies the 
Laplace transform of the time evolution of the centroid of the uncoupled (K = 0) ensemble 
to an infinitesimal displacement. Thus the stability of the coupled system is determined 
by properties of the collection of individual uncoupled elements. The formulation is valid 
for any type of dynamical behavior of the uncoupled elements. Thus ensembles whose 
members are periodic, chaotic, or a mixture of both can be treated. We discuss the analytic 
properties of M(s) and its numerical determination. We find that these are connected: 
analytic continuation of M(s) to the Im(s) axis is necessary for the application of the 
analysis, but, in the chaotic case (as discussed in Sees. |1V C| and [VT]) leads to numerical 
difficulties in obtaining M(s). We illustrate our theory by application to the Kuramoto 
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problem and by application to three different ensembles of globally coupled Lorenz systems. 
In particular, our Lorenz ensembles include a case where all of the uncoupled ensemble 
members are periodic with a pitchfork bifurcation of the uncoupled Lorenz equations within 
the parameter range of the ensemble, a case where all the ensemble members appear to be 
chaotic, and a case where the parameter range of the ensemble yields chaos with a window 
of periodic behavior. These numerical experiments illustrate the validity of our approach, 
as well as the practical limitations to numerical application. 
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Figures 



FIG. 1: Bifurcations. The horizontal line represents the incoherent state. Dashed (solid) lines 
represent unstable (stable) fixed points, and a system parameter governing the instability increases 
toward the right. 

FIG. 2: Illustraion showing how Mj is obtained for the Kuramoto example. 



FIG. 3: Bifurcation diagram for the Lorenz system in the parameter range r € [150, 165] (periodic 
ensemble) . 

FIG. 4: Periodic orbits from the Lorenz equation with (a) r = 165 and (b) r = 150. Black and 
grey in (b) denote the separate periodic orbits that are present after the pitchfork bifurcation. 

FIG. 5: The order parameter as a function of the coupling k for the periodic ensemble. Transitions 
are observed at (a) k* b = 0.95 ± 0.1 and (b) -\k* a \ = -0.70 ± 0.04. 

FIG. 6: Phase portraits showing the transition of the order parameter ((x(t))) for k slightly past 
(black) and slightly before (grey) the critical values at (a) — | /c * | and (b) k b . 

FIG. 7: Re[Mn(— iuj)] (black) and Re[Mn(— iu>)] (grey) vs. u, calculated using (a) the frequency 
response method (A = 0.05), and (b) the linear method described in Sec. IV Bj . Both methods 



yield very similar results overall. Predicted critical coupling values, calculated using the frequency 
response method and the values at the points indicated in (a), are — \k*\ = —0.72 ± 0.05 and 
k* b = 0.93 ±0.3. 

FIG. 8: The power spectrum of ((x^(t))) for k = 1.2, i.e. slightly larger than k^. The largest 
peak occurs at a frequency of 21.3 ± 0.1, in agreement with the predicted value 21.38 ± 0.15. 

FIG. 9: The (1,1) component of <(M(t)))» vs. t. 
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FIG. 10: The bifurcation diagram for the Lorenz system in the parameter range r <G [28, 52] (chaotic 
ensemble) . 

FIG. 11: The order parameter as a function of the coupling k for the chaotic ensemble. A subcritical 
transition is observed at — |fc*| = —5.56 ± 0.01. 

FIG. 12: Phase portrait of the transition of the order parameter (black). The central black oval is 
before the transition; afterwards, the order parameter shifts to one of the two lateral black dots. 
A single uncoupled Lorenz attractor for r = 52 is shown in the background (grey) for comparison. 

FIG. 13: Re[Mn(— iui)] (black) and Re[Mn(— iu)] (grey) vs. u. A = 2. For comparison, 
Re[Mn(— obtained with the linear (solid circles) and the impulse-response (open circles) meth- 
ods are included. 

FIG. 14: Graph of ((Mn (<)))* for the chaotic ensemble using the linear (grey) and the impulse 
response methods (black). 

FIG. 15: Graph of ln((x^(t))) vs. t showing the destabilization of the incoherent state for k 
slightly larger than k*. Since the transition is subcritical, only the initial growth rate is estimated 
as shown. 

FIG. 16: 7 vs. Sk for the chaotic ensemble near the transition. The slope predicted using the 
frequency response method is also shown (lines). 

FIG. 17: The bifurcation diagram for the Lorenz system in the parameter range r <G [167, 202] 
(mixed ensemble). 

FIG. 18: The order parameter as a function of the coupling k for the chaotic ensemble. Transitions 
are observed at -\k*\ = -0.68 ± 0.03 and jfe* = 1.75 ± 0.05. 

FIG. 19: Phase portraits showing the transition of the order parameter ((x(i))) for k slightly past 
(black) and slightly before (grey) the critical values at (a) — | k* | and (b) k£ . 

FIG. 20: Re[Mu(-iu)] (black) and Re[M u (-iu)] (grey) vs. u; A = 0.7 and N = 20,000. The 
predicted transition values are —\k*\ = —0.72 ± 0.09 and k^ = 1.64 ± 0.08. 
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FIG. 21: 7 vs. 5A; for the transition in the mixed ensemble near — \k*\. The slope predicted from 
the frequency response method is also shown (lines). 
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